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Some of our initial results show warning signs about extracting high-lying resonance energies 
using only single-hadron operators. 
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1. Introduction 



In a series of papers % |, |, [|| g], 



we have been striving to compute the finite-volume 



stationary-state energies of QCD using Markov-chain Monte Carlo integration of the QCD path 
integrals formulated on a space-time lattice. Such calculations are very challenging. The use 
of carefully designed quantum field operators is crucial for accurate determinations of low-lying 
energies. To study a particular state of interest, the energies of all states lying below that state 
must first be extracted, and as the pion gets lighter in lattice QCD simulations, more and more 
multi-hadron states lie below the masses of the excited resonances. The evaluation of correlations 
involving multi-hadron operators contains new challenges since not only must initial to final time 
quark propagation be included, but also final to final time quark propagation for a large number 
of times must be incorporated. The masses and widths of resonances must be deduced from the 
discrete spectrum of finite- volume stationary states for a range of box sizes. 

To compute the QCD stationary-state energies, the matrices Qj(t) of temporal correlations of 
sets of single-hadron and multi-hadron operators are estimated using the Monte Carlo method. For 
anNxN matrix, the N eigenvalues of C(t y 1 / 2 C(t)C(to)~ 1 / 2 tend to exp(-E k (t - 1 )) for large t 
and fixed to, where the decay rates are the N lowest-lying stationary-state energies that can be 
produced from the vacuum by the operators used. To compute the correlations involving isoscalar 
mesons and good multi-hadron operators, quark propagation from all spatial sites on one time slice 
to all spatial sites on another time slice are needed, and propagation from the sink time to the sink 
time are needed for a large number of sink times. A new method known as the stochastic LapH 
method^ has been introduced to make such computations accurate and practical in large volumes. 

2. The Stochastic LapH method 

Finding better ways to stochastically estimate slice-to-slice quark propagators is crucial to the 
success of our excited-state hadron spectrum project at lighter pion masses. We have developed 
and tested a new scheme which combines a new way of smearing the quark field with a new way 
of introducing noise and dilution projectors. The new quark-field smearing scheme, here called 
Laplacian Heaviside (LapH), has been described in Ref. [^] and is defined by 



where A is the three-dimensional covariant Laplacian in terms of the stout-smeared gauge field and 
o s is the smearing parameter. The gauge-covariant Laplacian operator is ideal for smearing the 
quark field since it is one of the simplest operators that locally averages the field in such a way that 
all relevant symmetry transformation properties of the original field are preserved. Let Va denote 
the unitary matrix whose columns are the eigenvectors of A, and let Aa denote a diagonal matrix 
whose elements are the eigenvalues of A such that A = Va Aa ■ The LapH smearing matrix is 
then given by 5? = Va & (of + Aa) v£. Let V s denote the matrix whose columns are in one-to-one 
correspondence with the eigenvectors associated with the N v lowest-lying eigenvalues of —A on 
each time slice. Then our LapH smearing matrix is well approximated by the Hermitian matrix 
5? = V s V/ . Evaluating the temporal correlations of our hadron operators requires combining Dirac 
matrix elements associated with various quark lines £2. Since we construct our hadron operators 




(2.1) 
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out of covariantly-displaced, smeared quark fields, each and every quark line involves the following 
product of matrices: 

g = D^yM~ 1 yDW\ (2.2) 

where DW is a gauge-covariant displacement of type i. 

An exact treatment of such a quark line is very costly and wasteful. Given our use of the Monte 
Carlo method to evaluate the path integrals over the gauge link variables, the statistical errors in 
our estimates of the hadron correlators are ultimately limited by the statistical fluctuations arising 
from the gauge-field sampling. Thus, we only need to estimate the quark lines to an accuracy 
comparable to the gauge noise from the Monte Carlo method. Such estimates can be obtained with 
far fewer inversions than required by an exact treatment of the quark lines. 

Random noise vectors tj whose expectations satisfy E{r\i) = and E^T)*-) = S;j are useful 
for stochastically estimating the inverse of a large matrix M as follows. Assume that for each of Nr 
noise vectors, we can solve the following linear system of equations: MX^ = r\ W for xW. Then 
X^> = M~ i r\( r \ and E{XiT\*j) = M^ 1 so that a Monte Carlo estimate of M^ 1 is given by M^ 1 ps 

lim^^oo 7^L^i^/ r ^j^*- Unfortunately, this equation usually produces stochastic estimates with 
variances which are much too large to be useful. Variance reduction is done by diluting the noise 
vectors. A given dilution scheme can be viewed as the application of a complete set of projection 
operators PH Define r]f ] = P^ffr, and further define X® as the solution oiM ik xf = t]l a \ then 
we have 

1 Nr 

The use of Z4 noise ensures zero variance in the diagonal elements £(t],-7]*). 

The effectiveness of the variance reduction depends on the projectors chosen. With LapH 
smearing, noise vectors p can be introduced only in the LapH subspace. The noise vectors p now 
have spin, time, and Laplacian eigenmode number as their indices. Color and space indices get 
replaced by Laplacian eigenmode number. Again, each component of p is a random Z4 variable so 
that E(p) = and E(pp J( ) = I^. Dilution projectors are now matrices in the LapH subspace. 
In the stochastic LapH method, a quark line on a gauge configuration is evaluated as follows: 

= D^yM- l v s v}D^\ 

= Zb d^^m- 1 v s p w v t 

= £ fo £)0')^M- 1 V s / J ( fe )£(pp t )/ J ( 6 ) t y/£)W t , 

= ZbE(p {i) yM- l V s P^p (D^VsPMp)^ . (2.4) 

For a noise vector labelled by index r, displaced-smeared-diluted quark source and quark sink 
vectors can be defined by 

p (r)[b](j) = D OV,pWpW, (2.5) 
pWMU) = D U)y M -i V s P w p {r \ (2.6) 
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and each quark line on a given gauge configuration can be estimated using 

^»+-tZd )m f# )m -, (2-7) 
N R r=l b 

where the subscripts u,v are compound indices combining space, time, color, and spin. 

Our dilution projectors are products of time dilution, spin dilution, and Laph eigenvector dilu- 
tion projectors. For each type (time, spin, Laph eigenvector) of dilution, we studied four different 
dilution schemes. Let N denote the dimension of the space of the dilution type of interest. For time 
dilution, N = N t is the number of time slices on the lattice. For spin dilution, N = 4 is the number 
of Dirac spin components. For Laph eigenvector dilution, N = N v is the number of eigenvectors 
retained. The four schemes we studied are defined below: 

P-j = §ij, a = 0, (no dilution) 

= S ai , a = 0, . . . , N - 1 (full dilution) 
pjf = Sij 5 8) [Kim a = 0,...,K-l, (block-*) 
Pif = Sij S a>im0 dK a = 0,. . . ,*- 1, (interlace-*) 

where i,j = 0,. . . ,N — I, and we assume N/K is an integer. We use a triplet (T, S, L) to specify a 
given dilution scheme, where "T" denotes time, "S" denotes spin, and "L" denotes Laph eigenvector 
dilution. The schemes are denoted by 1 for no dilution, F for full dilution, and B* and I* for 
block-* and interlace-*, respectively. For example, full time and spin dilution with interlace- 
8 Laph eigenvector dilution is denoted by (TF, SF, LI8). Introducing diluted noise in this way 
produces correlation functions with significantly reduced variances, yielding nearly an order of 
magnitude reduction in the statistical error over previous methods. The volume dependence of 
this new method was found to be very mild, allowing the method to be useful on large lattices. 
For all forward-time quark lines, we use dilution scheme (TF, SF, LI8), and for all same-sink-time 
quark lines, we use (TI16, SF, LI8). Interlacing in time makes the calculation of a large number of 
sink-to-sink diagrams possible with a feasibly small number of Dirac matrix inversions. The use 
of dilution projectors that interlace in time is perhaps the most important factor behind the success 
of the stochastic LapH method. 

Details on how the temporal correlations of hadron operators are evaluated are given in Ref. [gj. 
A very useful feature of the method is the fact that the hadron correlators completely factorize into a 
function associated with the sink time slice tp, and another function associated with the source time 
slice to- Summations over color, spin, and spatial sites at the source can be completely separated 
from the color, spin, and spatial summations at the sink. The stochastic LapH method leads to com- 
plete factorization of hadron sources and sinks in temporal correlations, which greatly simplifies 
the logistics of evaluating correlation matrices involving large numbers of operators. Implementing 
the Wick contractions of the quark lines is also straightforward. Contributions from different Wick 
orderings within a class of quark-line diagrams differ only by permutations of the noises at either 
the source or the sink. 

3. First results 

Our first results for the iso vector mass spectrum on a large 24 3 x 128 anisotropic lattice are 
shown in Fig. [j] The pion mass is about m % ~ 390 MeV here. These results are not finalized 
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Figure 1: Masses of the isovector mesons in terms of the nucleon mass using the stochastic LapH method 
with 170 gauge configurations on a 24 3 x 128 anisotropic lattice. The pion mass is about m K ~ 390 MeV. 
In the irrep labels, the letters with numerical subscripts refer to the point group (9/, irreps, the subscripts g 
and u refer to even and odd parity, respectively, and the superscripts ± refer to G-parity. Only single-hadron 
operators were used with dilution scheme (TF,SF,LI8). The shaded region indicates the threshold locations 
for multi-hadron energy levels. We emphasize that extractions of energies in the shaded regions could be 
complicated by "false plateaux" unless multi-hadron operators are used. 



since only single-hadron operators were used. The threshold locations for multi-hadron energy 
levels, assuming the interaction energies are small, are indicated by the shaded region. Extractions 
of energies in the shaded regions can only be considered reliable if multi-hadron operators are 
used, in addition to the single-hadron operators. The inclusion of the multi-hadron operators is in 
progress. 

Some initial results that incorporate multi-hadron operators are shown in Fig. [|. A 2 x 2 cor- 
relation matrix was evaluated involving a single-site p -meson operator and a total isospin I = \ TtTt 
operator in a P-wave with minimal relative momentum. The stochastic LapH method enables very 
accurate estimates of all elements (both diagonal and off-diagonal) of this correlation matrix such 
that diagonalization can be done. The effective masses associated with the diagonalized correlator 
are shown in Fig. |2[ The p operator dominates the lowest-lying level, while the %% operator domi- 
nates the first-excited state. Although the mixing of the operators is small, it is not negligible. This 
is certainly a warning about the dangers of extracting high-lying resonance energies using only 
single-hadron operators. 

Determining meson masses in the interesting scalar isoscalar sector will ultimately involve 
including a scalar glueball operator, so we began looking into the feasibility of such calculations. 
LapH quark-field smearing involves the covariant spatial Laplacian A. The eigenvalues of the 
Laplacian are invariant under rotations and gauge transformations so are appropriate for a scalar 
glueball operator. The lowest-lying eigenvalue was studied, as well as other functions of the eigen- 
values. We found that any combination of the low-lying eigenvalues worked equally well for 
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Figure 2: Effective masses corresponding to the diagonal elements of the rotated 2x2 correlator ma- 
trix involving a single-site p -meson operator and an / = 1 nn operator in a P-wave with minimal relative 
momentum. The p operator dominates the lowest-lying level (left), while the %% operator dominates the 
first-excited state (right). Although the mixing of the operators is small, it is not negligible. These results 
were obtained on 584 configurations of the 24 3 x 128 lattice with pion mass m K ~ 240 MeV. 
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Figure 3: Comparison of the effective mass associated with the correlator of the standard smeared plaquette 
glueball operator (left) with that of our new glueball operator Ga defined in the text (right). The similarity 
of the two results shows that the new glueball operator is just as useful as the familiar smeared plaquette for 
studying the scalar glueball. These effective masses do not reach a plateau at the glueball mass since various 
%% states and other multi-hadron states have smaller energies than the glueball mass. These results were 
obtained using 584 configs of the 24 3 ensemble for pion mass m K ~ 240 MeV. 



studying the scalar glueball. In particular, the operator defined by G\(t) = — Tr(@(a 4 2 + A)A(?)) 
was used. The effective mass associated with this operator is compared to that of the smeared pla- 
quette glueball operator in Fig. [|. The similarity of the results shows that Ga is just as useful as the 
familiar smeared plaquette for studying the scalar glueball. 

Effective masses corresponding to the diagonalized 4x4 correlator matrix involving a scalar 
isoscalar single-site quark-antiquark meson operator, the new glueball operator Ga, and two 1 = 
Tin operators in an S wave (one with zero relative momentum and the other with minimal nonzero 
relative momentum) are shown in Fig. |j. Mixing of these operators is sizeable. More %% operators 
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Figure 4: Effective masses corresponding to the diagonal elements of the rotated 4x4 correlator matrix 
involving a scalar isoscalar single-site quark-antiquark meson operator, the new glueball operator Ga, and 
two / = 7ZK operators in an S wave (one with zero relative momentum and the other with minimal nonzero 
relative momentum). Mixing of these operators is sizeable. More nn operators must be included to reliably 
extract the glueball mass. These results were obtained on 100 configurations of a 16 3 x 128 lattice with pion 
mass m K ~ 390 MeV. 



must be included to reliably extract the glueball mass. 



4. Conclusion 

A new method of stochastically estimating the low-lying effects of quark propagation was pro- 
posed which allows accurate determinations of temporal correlations of single-hadron and multi- 
hadron operators in lattice QCD. The method enables accurate treatment of hadron correlators 
involving quark propagation from all spatial sites on one time slice to all spatial sites on another 
time slice. Contributions involving quark lines originating at the sink time tp and terminating at the 
same sink time tp are easily handled, even for a large number of tp times. 

The effectiveness of the method can be traced to two of its key features: the use of noise dilu- 
tion projectors that interlace in time and the use of noise in the subspace defined by the Laplacian 
Heaviside quark-field smearing. Introducing noise in the LapH subspace results in greatly reduced 
variances in temporal correlations compared to methods that introduce noise on the entire lattice. 
Although the number of Laplacian eigenvectors needed to span the LapH subspace rises linearly 
with the spatial volume, we found that the number of inversions of the Dirac matrix needed for 
a target accuracy was remarkably insensitive to the lattice volume, once a sufficient number of 
dilution projectors were introduced. 

In addition to increased efficiency, the stochastic LapH method has other advantages. The 
method leads to complete factorization of hadron sources and sinks in temporal correlations, which 
greatly simplifies the logistics of evaluating correlation matrices involving large numbers of opera- 
tors. Implementing the Wick contractions of the quark lines is also straightforward. Contributions 
from different Wick orderings within a class of quark-line diagrams differ only by permutations of 
the noises at the source. 

The results presented here demonstrate that the stochastic LapH method is useful for accu- 
rately estimating all of the temporal correlations needed for a full study of the QCD stationary- 
state energy spectrum, which we are currently pursuing. A new glueball operator was tested and 
computing the mixing of this glueball operator with a quark-antiquark operator and multiple two- 
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pion operators was shown to be feasible. Some of our initial results showed warning signs about 
extracting high-lying resonance energies using only single-hadron operators. 
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